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Abstract 

We present results of fully non-linear calculations of decay of the inflaton 
interacting with another scalar field X. Combining numerical results for cos- 
mologically interesting range of resonance parameter, q < 10 6 , with analytical 
estimates, we extrapolate them to larger q. We find that scattering of X fluc- 
tuations off the Bose condensate is a very efficient mechanism limiting growth 
of X fluctuations. For a single-component X, the resulting variance, at large 
q, is much smaller than that obtained in the Hartree approximation. 
PACS numbers: 98.80.Cq, 05.70.Fh 
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In recent years, we have come to realize that the post-inflationary universe had probably 
been a much livelier place than was previously thought. In many models of inflation, the 
decay of the inflaton field is not a slow perturbative process but a rapid, explosive one. At 
the initial stage of this rapid process, called preheating [|| , fluctuations of Bose fields coupled 
to the inflaton grow exponentially fast, which can be thought of as "parametric resonance" 
@, and achieve large occupation numbers. At the second stage, called semiclassical ther- 
malization ||, the resonance smears out, and the fields reach a slowly evolving turbulent 
state with smooth power spectra ||,|J. 

The explosive growth of Bose fields leads to very large variances of these fields close 
to the end of the resonance stage. That could result in several important effects taking 
place shortly after the end of inflation. These include symmetry restoration, baryogenesis, 
and SUSY breaking []5| . To find out if these effects had indeed occurred, one needs a good 
estimate of the maximal size of Bose fluctuations. The semiclassical nature of processes 
involving states with large occupation numbers allows us to treat decay of the inflaton, 
and any Bose condensate in general, as a classical non-linear problem with random initial 
conditions for fluctuations fl3|. This classical problem can be analyzed numerically. 

In this Letter we report results of fully non-linear calculations for the most interesting 
case when the coupling of a massive inflaton <fi to some other scalar field X is relatively large. 
That means, more precisely, that the system is in the regime of wide parametric resonance 
|J, characterized by a large value of the resonance parameter q, q ^> 1. We have studied 
both expanding and static universes, to cover both post-inflationary dynamics and decays of 
possible other Bose condensates. Our objective was to obtain an estimate for the maximal 
size of X fluctuations, importance of which we emphasized above. 

Our results are as follows. We have found that scattering of X fluctuations off the Bose 
condensate of </>, which knocks inflaton quanta out of the condensate and into low-momentum 
modes, is very efficient in limiting the size of X fluctuations for large values of q, such as 
required [l],f§] to produce particles much heavier than the inflaton. This scattering process 
involves the condensate of zero-momentum inflatons and, for that reason, is especially en- 
hanced, cf. Refs. []3|||0]- Fluctuations of X can reach larger values for smaller values of q. 
The suppression of the maximal size of X fluctuations for large q significantly restricts the 
possibility of GUT baryogenesis after inflation, as well as the types of phase transitions that 
could take place after preheating. 

Our present results should be compared with those obtained in the Hartree approxi- 
mation. We find that for a single- component field X in flat space-time, the Hartree ap- 
proximation is inadequate for all q ^> 1. A similar conclusion was made in Ref. H] for the 
conformally invariant case of massless inflaton interacting with a massless field X, based on 
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our simulations of the fully non-linear problem for that case. The Hartree approximation, 
with its characteristic positive feedback of X on the inflaton decay |],U, may still apply 
when X has sufficiently many components; it remains to see if this can happen for realistic 
sizes of GUT multiplets. 

In the model with a massive inflaton (Model 1 of Ref. (|J), the full scalar potential is 
Vi(4>, X) = ^m 2 (j) 2 + ^g 2 (j) 2 X 2 + \M\X 2 . For comparison, we will sometimes present results 
also for the model with massless inflaton (Model 2), in which the potential is ^(0, X) = 
^A0 4 + |g 2 2 X 2 . Both fields (0 and X) have standard kinetic terms and are minimally 
coupled to gravity. We will often use rescaled variables: for Model 1, r = 777.77, where 77 is 
the conformal time; £ = mx; (p = 0a(r)/0(O); \ = Xa(r)/0(O). For Model 2, one should 
replace m with v^A0(O) in these rescalings. Here a(r) is the scale factor of the universe, 
defined so that a(0) = 1, and 0(0) is the value of the homogeneous inflaton field at the end 
of inflation. The resonance parameter q is q = g 2 2 (O)/4m 2 for Model 1, and q = c? 2 / 4A for 
Model 2. 

In order for resonance to fully develop in expanding universe, the resonance parameter 
q should exceed a certain minimum value, g m i n , which depends on the mass M x of X [|J. 
In Model 1, for Mx = 0, q m \ n ~ 10 4 ; for Mx = 10m, g min ~ 10 8 . We have simulated this 
model for q up to q = 10 4 in flat space-time and for a few q ranging from 10 4 to 10 6 in 
the expanding universe. We have developed an analytical approach, and we use analytical 
estimates to extrapolate our results to larger q, needed to produce heavier particles. 

The full non-linear equations of motion for (p and x, which follow from the action de- 
scribed above, are solved directly in configuration space with initial conditions corresponding 
to conformal vacuum at the end of inflation for all modes with non-zero momenta. The ini- 
tial conditions for x are given in Ref. J3J; those for S(p are obtained similarly, along the lines 
of Ref. . Classical fluctuations evolve from quantum ones, and, in cases we consider, their 
typical initial sizes are much smaller than the scale of non-linearity. Hence, initial evolution 
is linear with respect to fluctuations. During the initial linear stage, the equation of motion 
for Fourier components of x can be approximated as + w i( r )Xk = 0, where 

ulir) = m 2 x a 2 + k 2 -a/a + 4qip 2 {T) , (1) 

m x = Mx/m, and (fo is the zero-momentum mode of (p; <po(0) = 0. By virtue of our 
rescaling, ^0(0) = 1. 

Computations were done on 128 3 lattices, for a single-component X. For the case of 
expanding universe, we used m 2 = lO~ 12 Mp l . Energy non- conservation in flat space-time 
typically was less then 10~ 3 ; in expanding universe in linear regime our calculations closely 
reproduced calculations in the Hartree approximation, which have much better accuracy. 
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Let us first consider the case without expansion of the universe. In the formulas above, 
one substitutes a(r) = 1. The variances of the fields x an d y 9 in Model 1 at q — 2000 as func- 
tions of time are shown in Fig. |l[ The exponential growth of the variance (x 2 ) (the angular 
brackets denote averaging over space or, for space-independent quantities, equivalently, over 
realizations of random initial conditions) at early times is a parametric resonance, which in 
the present case ends at r m 40. This marks the end of the linear stage. 

At large q, fluctuations of X are produced at resonance stage only during short intervals 
of time near moments when <p passes through zero These are intervals in which the 

adiabatic (WKB) condition ci^/cc^ 3> 1 is broken for some k. Notice a series of spikes in (x 2 ) 
at the same moments of time. They are due to modulation of ujk by the oscillating tpo. Indeed, 
introduce analogs of occupation numbers via (x 2 ) = J d 3 kP x (k) oc / d^knk{T) / oj^ir) . 
Even at the resonance stage, change in n& during one oscillation is much smaller than 
variation of lo^. when instead of (x 2 { T )) we plot / d 3 krik(r), the spikes are replaced by 
relatively small steps at times when ipo = 0. 

We have monitored the power spectra, P(k), of <p and x m & U our integrations. The 
strongest resonant momentum of x is typically of order g 1 / 4 ; for some q, though, it can be 
close to zero. Development of resonance peaks for x is followed by appearance of peaks for ip 
due to rescattering, cf. Ref. ||. Later, rescattering leads to a turbulent state, characterized 
by smooth power spectra. 

Resonant growth stops when the Hartree correction to the mass of (p, due to (x 2 ) (at 
spike, since this is when fluctuations are produced) moves the system out of resonance. 
Because some modes grow faster than others, the width of the resonance peak in the power 
spectrum of x is much smaller than the full width (of order q 1 ^ 4 ) of the instability band. 
Nevertheless, we find, using analytical results of Ref. ||, that the peak width still scales 
as g 1//4 , up to a power of lng. We then estimate that resonance ends when (x 2 )s ~ q~ 3 ^ 2 . 
In Fig. ^ we have plotted results of our numerical integrations for the spike value (x 2 )s 
at the end of the resonance, as a function of q. At q > 100, these data are well fitted by 
(x 2 ) s q^ 3 ^ 2 ■ Figure [2] confirms that the termination of parametric resonance is a Hartree 
effect. 

Parametric resonance is followed by a plateau (unless we consider an exceptional q, for 
which the resonance peak was close to zero). There, the variances of fluctuations do not grow, 
but an important restructuring of the power spectrum of x takes place. The power spectrum 
of x changes from being dominated by a resonance peak at some non-zero momentum to 
being dominated by a peak near zero. The width of this new peak at k ~ is of order one. 
When this peak becomes strong enough, the growth of variances resumes. The resumed 
growth (at r > 50 in Fig. Q) is quite rapid (compared to subsequent slow evolution) and is 
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strongly affected by rescattering. This stage can be called semiclassical thermalization ||, 
or chaotization stage. Towards the end of it, the power spectra smoothen out; both power 
spectra are now dominated by momenta of order one. 

An important effect seen in Fig. [I] is the rapid growth of fluctuations of the field <p. 
Indeed, at late times, they are much larger than fluctuations of x- Fluctuations of p are 
produced by the scattering process in which x fluctuations knock ip out of the zero mode. If 
we neglected fluctuations of ip, i.e. the content of its non-zero modes, the Hartree approx- 
imation in our model would be exact. As we will see, at large q it can be far from being 
so. 

At the end of the chaotization stage (r « 80 in Fig. P, the variances in "valleys" between 
the spikes of (A 2 ) can be estimated directly from the classical equations of motion, solved 
in the limit \6(p\ <C <p~; here dip = ip — (fo, and (p is the amplitude of oscillations of ip . Using 
Green functions, we obtain, for Fourier components of 5ip and Xi 

M T) « "if £ M M%(r ~ r')]cp J d 3 k X ix k+p (2) 

Xk( r ) = *k + 8 <? [ dt ' F k(r, r') V o J rf 3 WpX k+P (3) 

where Q p = (p 2 + 1) 1//2 ; Xk ( T ) sorv es the Hartree equation, and Fk(r,r') is the retarded 
Green function for it. At times when x grows, the main contribution to the time integral 
in (H) comes from r' near r. Using the random phase approximation, we then obtain 
(Iv^pl 2 ) ~ q 2 <p 2 V~ 1 J rf 3 ^(|Xkl 2 ) (l^k+p| 2 )' wn ere V is the total spatial volume. The end of 
chaotization stage is the time, T c h, when the second term in r.h.s. of (|f) becomes of the 
order of the first. We find that that happens when, in terms of the physical fields, 

^ ~ 3/w 0) S i ((50) 2 )~0 2 (O)/g. (4) 

where is the amplitude of O - I n particular, (A 2 ) in "valleys" after rescattering, Eq. (|J), 
is of the order of (A 2 ) s at the end of resonance. This is indeed seen in Fig. [IJ as well as in 
all our other integrations in flat space-time. Note that when variances of fields reach values 
(|j), the analogs of occupation numbers, introduced as before, reach values of order 1/g 2 for 
both X and 0. 

After chaotization, the maxima and minima of (x 2 ) evolve slowly, but the variance of ip 
continues to grow rapidly for a while, see Fig. [I]. We interpret this as follows. According to 
Eq. (H), fluctuations of <p are driven by x- A periodic (or close to periodic) x causes growth 
of 5(p via a usual, i.e. not parametric, resonance. This growth of 5<p will stop only when the 
approximation leading to (0) breaks down, that is 8(f) ~ <fi- 
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Let us now turn to Model 1 in expanding universe, see Figs. We used 0(0) = 0.28Mpi 
For massless A, the evolution of the scale factor was determined self-consistently, in- 
cluding the influence of produced fluctuations in the Einstein equations; for massive A, the 
universe was assumed matter dominated. 

In expanding universe, particle creation acquires a qualitatively new feature |§] : because 
of the time-dependence of (p, the resonance peak scans the entire instability band, see Fig. |j. 
As our numerical integrations confirm (see also |4j]), in order for production of fluctuations 
to be efficient, variation of the frequency of x, Eq. (P, need not be periodic. What is 
required is that every once in a while the adiabatic condition breaks down. In this situation, 
"non-adiabatic amplification" seems to be a better term for stimulated particle creation than 
"parametric resonance" . 

The variances of fields at the end of chaotization can be estimated in the same way as 
before, and we obtain again Eq. (|J). Note that now 4>{r c h) <C 0(0), due to the redshift of 
the field. 

Evolution of (X 2 ) after chaotization can be thought of as a result of balance between 
creation and destruction of X fluctuations. These processes are much faster than the ex- 
pansion of the universe at this stage. We may then assume that the system looses memory 
of initial conditions and evolves as it would in flat space-time but with the effective, slowly 
changing value of the resonance parameter, q c fi(r) = g0 2 (r)/0 2 (O). This assumption allows 
us to extend the estimate @ for (X 2 ) v to times after chaotization by simply using the cur- 
rent value of in it. In general, the assumption allows us to relate the scaling of variances 
with to their scaling with q: (A 2 )/0 2 oc 9^ a (r) oc q~ a 0~ 2a (r) . 

Time dependence of variances is shown in Fig. |3] for the case q = 10 6 and m x = 2. We 
can construct an analog of particle density, n, as n/m = [4g(0Q + (50) 2 )/0 2 (O) + m 2 ] 1//2 (A 2 ). 
As long as the term with O is the leading term in the square brackets, n oc (A 2 )„0, 
and according to our scaling argument [(X 2 ) v scales with a = 3/2, cf. Eq. (|)], n after 
chaotization should be time-independent. Our data confirms that. Thus, fluctuations keep 
being produced at the same rate as they are diluted by the expansion. This can be called 
negative feedback amplification. A striking feature of the data is that, after chaotization, 
(A 2 ) in spikes becomes essentially time- independent. According to our scaling argument, 
this means that (A 2 ) s scales as 1/q. In Fig. ^|, we show by stars the values of (A 2 ) in 
expanding universe, taken at times when O decays (more precisely, when 2 « (<50 2 )) and 
the "spike" and "valley" values of (A 2 ) coalesce. (These values are within a factor of few 
from the maximal values (A 2 ) achieves in the expanding universe; the exact value is model 
dependent and it was largest, equal to 2.6 , for Mx = 0, q = 10 6 .) The scaling of these 
values with q is well fitted by 1/q. 
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We can now use the 1/q scaling of (X 2 ) s to extrapolate our data to larger values of q. For 
example, for q = 10 8 , required to produce scalar leptoquarks with Mx = 10m, we obtain 
(X 2 ) s ~ 10~ 10 Mp l . Because (X 2 ) s stays on this level for a while, despite the expansion, 
the time-integrated conversion of inflatons into X fluctuations in this two-field model is in 
fact quite efficient. For definite predictions for baryon asymmetry generated in decays of 
leptoquarks, however, one has to include other fields, which can have smaller q and thus 
provide faster alternative channels of inflaton decay. 

If instead of a single-component X we consider X with N real components, our estimate 
for (X 2 ) v at the end of chaotization becomes larger by a factor of ^/N. For realistic GUT 
values of N, this can increase the maximal value of (X 2 ) by an order of magnitude. In 
comparison, in the Hartree approximation, the maximal value of (X 2 )/(j) 2 (0) is of order 
<T 1/2 0/0(O) for q = 10 s and M x = 10m we get (X 2 ) max ~ lO^M^. 

The relatively quick complete exponential decay of (fi , seen in Fig. [3|, is a distinctive 
feature of Model 1, as opposed to the conformally invariant Model 2. In the latter case, 
instead of (§]), we obtain, at the end of chaotization, ((£0) 2 ) ~ 4> 2 /q, (X 2 ) v ~ <fi 2 /q 3 ^ 2 , and 
in the subsequent evolution (X 2 ) redshifts together with (j) 2 . 
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FIG. 1. Variances of the fields X (solid curve) and 4> (dotted curve) in Model 1 in flat space-time. 
The filled square marks the spike value of (X 2 ) at the end of resonance. 
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FIG. 2. Filled squares and crosses are the spike values, (X 2 ) s , at the end of resonance obtained 
in fully non-linear simulations of Model 1 (flat space-time) and Model 2; empty boxes are the spike 
values at the first plateau in the Hartree approximation. Stars correspond to (X 2 ) at the moment 
when zero mode decayed in Model 1 in expanding universe. 
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